Method for characterising particles in supension from frequency domain photon migration measurements

ABSTRACT

Methods are provided for measuring isotropic scattering coefficients of suspensions using multiply scattered radiation that is modulated in amplitude at selected modulation frequencies. The radiation may be light. Quantities describing diffusion of the multiply scattered radiation are preferably measured at a plurality of distances between source and receiver and a plurality of frequencies. Linear regression techniques are provided for maximizing accuracy of the scattering data at a selected wavelength of a radiation. Methods are provided for inversing an integral equation so as to determine a calculated value of scattering coefficient. Parameters are varied to minimize the difference between the calculated and measured scattering coefficients and thereby to determine volume fraction, particle size distribution and interparticle force between the particles in a suspension. By incorporating a first principles model to account for interparticle force, the measurements can be used to determine a parameter governing interparticle forces in a suspension. The suspension may be in a liquid or a gas.

The U.S. government has a paid-up license in this invention and theright in limited circumstances to require the patent owner to licenseothers on reasonable terms as provided for by the terms of Grant No.CTS-9876583, awarded by the National Science Foundation.

This invention relates to measurements of particle characteristics usingmultiply scattered radiation. More particularly, methods for obtainingmore precise Frequency Domain Photon Migration (FDPM) data and obtainingfrom the data particle size distributions, volume fractions, andinterparticle force characteristics for interacting particles andpolydisperse or multi-modal size distributions are provided.

On-line measurements of particle size distribution and volume fractionof colloidal suspensions are needed in various processes of chemical andother industries. The measurements may be used, for example, toevaluate, optimize, and control manufacturing processes or to controlproduct quality.

On-line particle-sizing methods such as turbidity, angular staticscattering and dynamic light scattering, based on light scattering froma single particle, are often used. These methods require considerabledilution of most commercial suspensions to avoid the effects of multiplescattering and particle interaction. Dilution may alter the sizedistribution of the particles. Otherwise used without dilution,empirical calibration of the instruments with like-suspensions ofparticles to be sized is necessary. Methods have been proposed tosuppress the effects of multiple scattering, such as fiber optic dynamiclight scattering (FODLS) and modified static light scattering, but thesemethods are not well suited for particle sizing of multiple scattering,dense colloidal suspensions.

Recently, new optical particle sizing methods, such as Diffusing WaveSpectroscopy (DWS) and diffuse reflectance and/or transmittancespectroscopy have been used for dense colloidal suspensions. To provideaccurate sizing information these methods require accounting forparticle interactions at high volume fractions.

In recent years, alternatives to optical sizing methods have beendeveloped. These include acoustic and electroacoustic spectroscopy.Similarly to light scattering methods, the acoustic and electroacoustictheories for dilute systems are relatively complete, but the theoriesfor polydisperse concentrated systems are far from complete. Whileseveral models based on first principles have been developed to accountfor particle interactions in concentrated suspensions, few of the modelshave been successfully used in practical applications. Also, theseacoustic methods require knowledge of the physical properties of eachphase and electroacoustic spectroscopy can only be used for chargedparticles. Better methods for characterizing dense suspensions areneeded. Further background information on these methods can be found inthe article “Approach for Particle Sizing in Dense PolydisperseColloidal Suspension Using Multiple Scattered Light,” Langmuir 2001, 17,6142-6147, which is hereby incorporated by reference herein.

Frequency Domain Photon Migration (FDPM) has been proposed for particlemeasurements (U.S. Pat. No. 5,818,583, for example). Since thistechnique measures time-dependent propagation characteristics ratherthan the amount of light detected, it has the major advantage of notrequiring calibration. In addition, FDPM measurements allowdetermination of absorption and scattering properties independently.Therefore, FDPM allows characterizations of colloidal suspension byscattering measurements that are not biased by changes in lightabsorption or color of suspending fluids. Because FDPM depends onmultiply scattered light, it is restricted to particulate suspensionsthat are concentrated enough to cause multiple scattering events asradiation is transmitted through the suspension. Such suspensions aretypically opaque. Using the FDPM method, both size distribution andvolume fraction can be directly obtained by inversion algorithms, wherethe only physical parameter required a priori is the relative refractiveindex between the particulate and continuous phases. Since no dilutionis necessary to avoid multiple scattering, FDPM is suitable for onlinemonitoring of particle size distribution and volume fraction in manycommercial processes that include particles or droplets. FDPM has beensuccessfully used for recovery of particle size distribution and volumefraction in opaque, multiple scattering suspensions of polystyrene andof titanium dioxide, for example (S. M. Richter et al., “Particle SizingUsing Frequency Domain Photon Migration,” Part. Part. Syst. Charact.,15, 9 ,1998).

Using visible wavelengths, particles from 50 nm to greater than 1 micronin diameter may be characterized in multiple scattering, non-interactingsuspensions. However, at volume fractions larger than about 1-5%,corrections for particle interaction become increasingly necessary.Particle interactions can arise from volume exclusion effects as well asfrom particle interactions, which arise owing to forces that act betweenparticles. As a results of these interactions, the particle suspensionbecomes ordered and structured. This order and structure impactsscattering and must not only be used to perform sizing, but when sizinginformation is available, can be used to determine the nature ofparticle interactions, (S. M. Richter et al, “Characterization ofConcentrated Colloidal Suspensions Using Time Dependent Photon MigrationMeasurements,” Colloid Surf. A, 172, 163, 2000).

Determination of particle characteristics from FDPM data requiressolving mathematically an “ill-posed” inverse problem. Error in the FDPMdata used to solve the inversion equations has a large effect onparticle size distribution, particle size, and other calculated results.Therefore, there is a need for apparatus and method for determining theFDPM (scattering) data to a high degree of accuracy. Further, a methodfor indicating the quality of the data is needed, such as with a qualityparameter that indicates whether the measurement is precise enough foruse in an inverse solution. In application Ser. No. 09/297,895, which ishereby incorporated by reference, measurements at multiple modulationfrequencies and two locations were disclosed. Since results arenon-linear at multiple frequencies, a non-linear regression method ofanalysis was necessary. A linear method of analysis and method fortreating results of measurements from multiple measurements for maximumaccuracy are needed.

In previous work on FDPM, particle size distributions were assumed to beeither normal, log-normal or Weibull distributions. A method fordetermining an unknown and arbitrary particle size distribution isneeded for polydisperse and multi-modal suspensions, at non-interactingor interacting (high) volume fractions.

Many suspensions of commercial interest contain volume concentrations ofthe disperse phase far above the 1-5 per cent by volume range, whereparticles can be considered non-interacting. In addition to the methodsfor determining volume fraction and particle size distribution whenparticle interaction must be considered, methods for indicating anyattractive or repulsive forces between the interacting particles areneeded, since such forces are important in determining stability,rheological properties, and other properties of the high-volume-fractionsuspensions. These interparticle forces, along with the volume exclusioneffects, cause additional order within the suspension. Therefore,methods for interpreting interparticle forces as result of the order areneeded. In applying FDPM to a particular sample containing particles,various models may be used to interpret a measured scatteringcoefficient in terms of particle concentration, particle sizedistribution and inter-particle forces. In the method disclosed in Ser.No. 09/297,895, interparticle interactions arising from volume exclusionwere interpreted assuming a monodisperse sample. Methods for consideringmore general particle size distributions are needed. A device thatincludes a processor to execute computer programs to provide results interms of the various models is needed.

FIG. 1(a) depicts a radiating wave of varying intensity propagatingthrough a scattering medium and FIG. 1(b) illustrates source anddetected waves along with measurements used in FDPM measurements.

FIG. 2(a) shows an isometric view and FIG. 2(b) a plan view of oneembodiment of apparatus for FDPM measurements.

FIG. 3 shows a schematic of one embodiment of instrumentation that maybe used in FDPM measurements.

FIG. 4 illustrates a flowchart showing one embodiment of a method fordetermining scattering coefficients from FDPM measurements.

FIG. 5 shows an illustration of the values of the criteria parameter, Q,at different frequencies and distances from source to detector.

FIG. 6 shows an “L-curve” used for parameter-choosing method forregularization of the solution of the inversion problem for particlesize distribution.

FIG. 7 illustrates a flowchart showing one embodiment of a method fordetermining particle size distribution and volume fraction of asuspension of non-interacting particles from scattering coefficients.

FIG. 8 illustrates a flowchart showing one embodiment of a method forsolving the non-linear inverse problem for any arbitrary model ofinteraction between particles, whether purely volume exclusion orincluding attractive and repulsive forces.

FIG. 9(a) and FIG. 9(b) show a flowchart of the linear approachassociated with the LMA approximation for determining particle sizedistribution, volume fraction and interaction parameter for a suspensionof particles which interact through volume exclusion.

FIG. 10(a) and FIG. 10(b) show a flowchart of the linear approachassociated with the DA approximation for determining particle sizedistribution, volume fraction and interaction parameter for a suspensionof particles which interact through volume exclusion.

FIG. 11 shows a comparison of different structure factor models withFDPM data at a wavelength of 687 nm.

FIG. 12 shows the change of FDPM-measured scattering coefficients withvolume fractions at two wavelengths.

FIG. 13 shows the changes of precision of FDPM-measured scatteringcoefficients with volume fraction at two wavelengths.

FIG. 14 shows the change of angle-integrated structure factors withvolume fractions at two wavelengths for latex samples.

FIG. 15 shows the change of scattering coefficient with volume fractionfor a polydisperse latex sample at a wavelength of 687 nm, showingmeasurements and calculated results with and without consideringparticle interaction.

FIG. 16 shows the effect of salt content in a sample on measuredscattering coefficients.

FIG. 17 shows the computation of an effective charge on particles ascationic surfactant is adsorbed onto negatively charged particles insuspension.

Method for measuring isotropic scattering coefficient for multiplyscattered radiation passing through a suspension is provided. In oneembodiment, the radiation is modulated at a selected frequency andmeasurements selected from among AC attenuation, DC attenuation andphase change of the radiation are made at different distances apart of asource and a detector of the radiation. In another embodiment, themeasurements are repeated at different modulation frequencies. Themeasurements may be used in linear equations to calculate absorption andscattering coefficient of the radiation. Simultaneous regression ofcombinations of the characteristics of the radiation is provided. A“criteria parameter” is provided for determining the experimentalquality of the measurements.

Method for determining volume fraction, particle size distribution andinterparticle force is provided. Differences between experimental valuesof scattering coefficient at selected wavelengths of radiation andcalculated values of scattering coefficient at the same selectedwavelengths are minimized using a least squares minimization. Theinversing of an integral equation describing calculated values ofscattering coefficient may be stabilized using Tikhonov regularizationand a non-negativity constraint. Particle size distribution and volumefraction are separated by a normalized distribution constraint.Determination of particle size distribution is not limited to a prioriassumptions. Monodisperse and polydisperse suspensions, either single-or bi-modal may be investigated using the methods disclosed herein.

The Hard Sphere model of particle interaction is applied, along withmodifications of this model. Models accounting for interparticle forcesare developed, including the inclusion of a first principles model ofparticle interaction, which makes possible investigation ofinterparticle forces that are important in determining rheology,stability and other properties of colloidal suspensions.

Apparatus is provided for measuring scattering coefficients thatincludes a source or multiple sources and a receiver or multiplereceivers at multiple distances apart in a sample of suspension. In oneembodiment, a source or detector is mechanically moved to allowmeasurements are multiple distances apart of the source and detector. Asystem for executing a computer program to perform the calculations isalso provided.

DESCRIPTION OF PREFERRED EMBODIMENTS

FIG. 1(a) depicts variations of light intensity in a scattering mediumbetween a modulated source and a detector. Frequency Domain PhotonMigration (FDPM) measurements involve the time-dependent propagationcharacteristics of multiply scattered light having a selected wavelengthin the scattering medium. U.S. Pat. No 5,818,583 and PCT application WO98/20323, which are incorporated by reference herein, disclose apparatusand methods for making and interpreting FDPM measurements for thepurposes of characterization of particles. The technique depends onintroducing intensity-modulated light into a multiply scattering mediumat a point source and detecting transmitted light at one or morediscrete points removed from the source. FIG. 1(b) illustratesquantities that are measured in obtaining FDPM data. Phase shift (PS)between-the source and detected signal, the amplitude of the alternatingsignal of the modulation at the source (ACs) and at the detector (ACd)and the time-average signal at the source (DCs) and at the detector(DC_(d)) may be measured. Measurements of these quantities are used asdescribed below to calculate properties of a scattering sample betweenthe source and detected signal.

FIG. 2 illustrates one embodiment of apparatus 10 that may be used forFDPM measurements. Experimental conditions preferably provide in effectan infinite scattering medium around a source of radiation, althoughother geometries are possible. This may be achieved by placing a fiberoptic source of radiation and fiber optic receiver probes in a samplecontainer. The container may have a volume of 100 mL, for example.Alternatively, one or more sources of radiation may be placed directlyin the sample and one or more detectors may also be placed directly inthe sample. The source may be a laser, an LED, or other source that canbe modulated in intensity in the frequency range where measurements areto be made. Preferably, a laser diode is used. The source may beinterchangeable or, adjustable, or be comprised of an array to produce aplurality of selected wavelengths of light. The source and one or morereceiver fibers may be placed at fixed distances apart or the relativeseparations between the source and detector fibers may be adjustedmechanically. Distances between source and detector fibers should be atleast ten isotropic scattering mean free paths, to satisfy conditions ofmultiple light scattering. Preferably, distances between the source andreceiver fibers are measured to an accuracy of at least about 0.1 percent.

Referring to FIG. 2, radiation source 14 provides modulated radiationinto source fiber 16, which ends in sample container 12. First receiverfiber 18(a), preferably ending in the same plane as the first fiber andat a spaced-apart distance, r_(o), in that plane is used to detect asignal propagating through the sample and to convey the signal tophotodetector 20. Second receiver optical fiber 18(b) may be placed at asecond spaced-apart distance, r₁, from source fiber 16. Additionalreceiver optical fibers 18(c, d, . . . ) may be placed at selectedspaced-apart distances from fiber 16, preferably in a pattern aroundfiber 16 so as to avoid placement of any fiber in the direct path ofradiation to another fiber. Alternatively, fibers may be moved from afirst distance, such as r_(o), to a selected number of distances, r,from the source. In another embodiment, source fiber 16 may be moved tomultiple locations for measurements at different source-receiverdistances while receiver fibers remain stationary.

Referring to FIG. 3, FDPM measurements may be performed with theelectronic apparatus illustrated. FDPM instrumentation is made up ofthree major components: (a) the light source, (b) the light delivery anddetection apparatus, and (c) the data acquisition and analysis system. A687 nm or 828 nm wavelength laser diode (Thorlab, Inc., Newton, N.J.)may be used as a light source, for example. Other lasers, LEDs or othersources of radiation may be used. The output of the laser diode may besinusoidally modulated at the frequencies of interest by a frequencysynthesizer (Marconi Instruments signal Generator 2022A, Mountain View,Calif.). Common modulation frequencies are in the range from 10 MHz to100 MHz. The modulated light from the laser source may be split by aglass beam-splitter so that approximately 10% of the light is collectedby a reference Photo-Multiplier Tube (PMT) (Model H6573, Hamamatsu,Japan) via one 1000 μm optical fiber (FT-1.0-EMT, Thorlab, Inc., Newton,N.J.). The remaining transmitted light may be launched into thescattering solution (sample) using a second 1000 μm source optical fiberpositioned vertically in the sample container. The scattered light fromthe sample may be detected by a third vertically positioned 1000 μmoptical fiber and directed to the sample photo-detector (PMT). Both PMTsmay be gain-modulated by an amplified RF signal from another frequencysynthesizer (Marconi Instruments Signal Generator 2022A, Mountain View,Calif.) which is phase-locked to the first. A standard heterodynetechnique may be employed to accurately process high frequency signals.Labview 5.0 software (National Instruments, Austin, Tex.) on a standardPC may be used as the platform for experimental data acquisition andanalysis. More detailed description of the apparatus may be found in thepaper “Precise Analysis of Frequency Domain Photon MigrationMeasurements for Characterization of Concentrated ColloidalSuspensions,” by Sun, Z., Huang, Y. and E. M. Sevick-Muraca, Rev. Sci.Instrum., 73(2): 383-393, 2002.

In concentrated particulate suspensions, the propagation of light can bemodeled as a diffusion process, since photons are multiply scattered.The governing equation for the propagation of light through a multiplescattering medium in the diffusion approximation can be written as:$\begin{matrix}{{\frac{\partial{U\left( {r,t} \right)}}{\partial t} - {{vD}{\nabla^{2}{U\left( {r,t} \right)}}} + {v\quad\mu_{a}{U\left( {r,t} \right)}}} = {q_{0}\left( {r,t} \right)}} & (1)\end{matrix}$where, U is the total photon density; q₀ is the isotropic source term; vis the speed of light in the medium; and D is the diffusion coefficient,which is defined as $\begin{matrix}{D = \frac{1}{3\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}} & (2)\end{matrix}$Herein, μ_(a) is the absorption coefficient and μ_(s)′ is the isotropicscattering coefficient that arises from multiple scatteringμ_(s)′(λ)=(1−g)μ_(s)(λ)   (3)where μ_(s) is the scattering coefficient, and g is the average of thecosine of the scattering angle, or the scattering anisotropy. (A list ofdefinitions of mathematical symbols used herein is provided at the endof this written description.)

In the case of a sinusoidally intensity-modulated point source of lightand a uniform infinite medium boundary condition, Equation (1) (alsocalled the photon diffusion equation) can be solved to yield analyticalexpressions for three experimentally determined quantities: (1) thesteady-state photon density or the time invariant average intensity,termed the DC component; (2) the amplitude of the photon-densityoscillation, termed the AC component; and (3) the phase shift of thephoton-density wave, termed PS. FIG. 1(b) illustrates those quantities.The solution to Eq. (1), as further described in the paper “Approach forParticle Sizing in Dense Polydisperse Colloidal Suspension UsingMultiple Scattered Light,” Langmuir 2001, 17, 6142-6147, which providesan expression for the isotropic scattering coefficient for apolydisperse interacting colloidal suspension, which is: $\begin{matrix}{{\mu_{s}^{\prime}(\lambda)} = {\int_{0}^{\infty}{{\frac{12\phi\quad{f(x)}}{k^{2}x^{3}}\left\lbrack {\int_{0}^{\pi}{{F\left( {n,x,\lambda,\theta} \right)}{S\left( {q,\phi} \right)}\sin\quad{\theta\left( {1 - {\cos\quad\theta}} \right)}\quad{\mathbb{d}\theta}}} \right\rbrack}\quad{\mathbb{d}x}}}} & (4)\end{matrix}$where F(n,x,λ,θ) is the form factor for a particle of diameter, x, withrelative refractive index of particle to medium, n, at wavelength, λ,and θ is the scattering angle. The term φ represents the volume fractionof particles in the suspension and f(x) represents Particle SizeDistribution (PSD); k is given by 2 πm/λ, where m is the refractiveindex of medium; q is the magnitude of the wave vector; q=2 k sin(θ/2).In equation (4), the form factor, F(n,x,λ,θ), can be calculated byclassical Mie scattering theory of single particle scattering, which iswell known in the art. The structure factor, S(q,φ), is a direct measureof the local ordering of colloidal particles, and the value of S(q,φ) isequal to unity in the absence of particle interactions (e.g. in a dilutesuspension).

For polydisperse systems, since particle interactions occur betweencolloids not only with the same particle size but also with differentparticle sizes, the partial structure factor, S_(i,j), must beconsidered, where i and j represent particles with different sizes x_(i)and x_(j). Equation (4) becomes: $\begin{matrix}{{\mu_{x}^{\prime}(\lambda)} = {\int_{0}^{\pi}{\frac{12\phi}{k^{2}}{\int_{0}^{\infty}{\frac{f\left( x_{i} \right)}{x_{i}^{3}}{\int_{0}^{\infty}{\frac{f\left( x_{j} \right)}{x_{j}^{3}}{{F_{i,j}\left( {n,x_{i},x_{j},\lambda,\theta} \right)} \cdot {S_{i,j}\left( {n,x_{i},x_{j},\lambda,\theta} \right)}}\sin\quad{\theta\left( {1 - {\cos\quad\theta}} \right)}\quad{\mathbb{d}x_{j}}\quad{\mathbb{d}x_{i}}\quad{\mathbb{d}\theta}}}}}}}} & (5)\end{matrix}$where F_(i,j) is the binary form factor between the particles withdifferent sizes x_(i) and x_(j).

If particle sizes are the same (i.e., i=j), the form factor is simplycalculated as a monodisperse system:F _(i,i)=(f _(1,i) f _(1,i) *+f _(2,i) f _(2,i)*)   (6)where f_(1,i) and f_(2,i) are the scattering amplitudes into twoorthogonal polarization states arising from a particle with size x_(i),which can be calculated by Mie scattering theory. f_(1,i)* and f_(2,i)*are conjugates of f_(1,i) and f_(2,i) respectively. However, if particlesizes are different (i.e., i≠j), the binary form factor with differentparticle sizes is calculated by:F _(i,j) =Re(f _(1,i) f _(1,j) *+f _(2,i) f _(2,j)*)   (7)

Referring to Eq. (5), since the scattering coefficient is a function ofwavelength, λ, measurements of scattering coefficient at M differentwavelengths provide a means for establishing a system of M equationsthat can be used to establish values for M parameters used to describescattering in a suspension. If measurements are made at only onewavelength, the number of unknowns must be reduced to one by assumptionsor other known conditions in the suspension. For example, if the volumeconcentration, φ, is low enough to neglect interparticle forces, thestructure factor, S, can be assumed to have a value of one. There arestill two unknowns, φ and f(x). If volume fraction, φ, is known and itis known that the sample is monodisperse, then the single particlediameter can be determined with measurements at only one wavelength, λ.If the sample of known volume fraction is polydisperse but assumed tohave a particle size distribution described as Gaussian or log-normal,which is described by two parameters, measurements at two wavelengthswill be necessary and sufficient to determine the two parameters used todescribe those distributions. If the particle size distribution is ofunknown shape, even possibly bi-modal, more complex size distributionmodels are necessary, as will be described below and more measurementsat different wavelengths will be necessary to describe the particle sizedistribution. If the volume fraction is not known, then measurements arerequired at an additional wavelength to determine the additional unknownparameter, φ, in Eq. (5).

It should thus be clear that Eq. (5) can be used for two general typesof investigations of suspensions: particle sizing and assessment ofparticle interactions. Methods are disclosed herein for particle sizingmeasurements when the form of the particle size distribution is known,which may be uni-modal or multi-modal, with or without particleinteraction. Methods are also disclosed for determining an unknownparticle size distribution, which may be uni-modal or multi-modal, withor without considering inter-particle interactions. These methodsinclude the assumption that the inter-particle interactions can beadequately treated using a “hard sphere” model or modifications of thatmodel to approximate the volume exclusion effects that predominantlyimpact the structure. To use Eq. (5) for assessment of particleinteractions, unknown or known forms of particle size distribution maybe used. An iteration parameter may be used to indicate repulsive orattractive interparticle forces on an empirical basis and without aspecified particle interaction model.

To determine the most accurate particle sizing data and particleinteraction assessment, accurate measurements of isotropic scatteringcoefficients at each wavelength must be available. A preferred methodfor determining isotropic scattering coefficient, μ_(s)′, is as follows:(1) using a source emitting a selected wavelength of light in the rangefrom about 100 nm to about 2000 nm, measure the phase shift and eitherthe AC or DC value of radiation intensity at a series of differingdistances from the source of light, preferably using a selectedmodulation frequency in the range from about 10 MHz to about 300 MHz.The number of differing distances is preferably at least three and ismore preferably in the range from about 5 to 10. Measurements atdifferent distances, called Multi-Distance (MD), are preferably repeatedat a series of modulation frequencies, called Multi-Frequency (MF).Preferably at least about three modulation frequencies are used, such as30, 60 and 90 MHz.

The following equations can conveniently be used to perform calculationsfrom the measurements: $\begin{matrix}{{\ln\left( {\frac{r}{r_{0}}{DC}_{rel}} \right)} = {- {\left( {r - r_{0}} \right)\left\lbrack {3{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}} \right\rbrack}^{1/2}}} & \left( {8a} \right) \\{{\ln\left( {\frac{r}{r_{0}}{AC}_{rel}} \right)} = {{- \left( {r - r_{0}} \right)}{\sqrt{3{{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}/2}}\left\lbrack {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} + 1} \right\rbrack}^{1/2}}} & \left( {8b} \right) \\{{\ln\left( {Mod}_{rel} \right)} = {{- \left( {r - r_{0}} \right)}{\sqrt{3{{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}/2}}\left\lbrack {\sqrt{2} - \left( {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} + 1} \right)^{1/2}} \right\rbrack}}} & \left( {8c} \right) \\{\left( {PS}_{rel} \right) = {\left( {r - r_{0}} \right){\sqrt{3{{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}/2}}\left\lbrack {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} - 1} \right\rbrack}^{1/2}}} & \left( {8d} \right)\end{matrix}$where AC_(rel) is the ratio of AC measured at r to that measured atr_(o), PS_(rel) is the difference in measured phase at r and r_(o), andDC_(rel) is the ratio of DC measured at r to that measured at r_(o).

Using equations 8a-8d, the optical absorption and scattering propertiesof the medium, μ_(a) and μ_(s)′, can be estimated from the relativemeasurement data at two or more source-detector distances. In addition,the following expression can be derived from the above equations$\begin{matrix}{\frac{{\ln^{2}\left( {{rAC}_{rel}/r_{0}} \right)} - {PS}_{rel}^{2}}{\ln^{2}\left( {{rDC}_{rel}/r_{0}} \right)} = {Q = 1}} & (9)\end{matrix}$Eq. 9 can be used to check whether FDPM measurements are consistent withthe assumptions of the photon diffusion model. We term Q the “criteriaparameter”. Since equation (9) is independent of the optical proprietiesμ_(a) and μ_(s)′, the value of Q can be used to optimize measurementconditions even before extracting the isotropic scattering coefficientby different data analysis approaches. So, it provides a usefulcriterion to assess the accuracy and precision of FDPM measurementsunder different experimental conditions.

A flowchart describing a preferred procedure for determining μ_(s)′ andμ_(a) from FDPM measurements as a selected wavelength of light is shownin FIG. 4. In step 402, a first modulation frequency is chosen, which isnormally in the range of 10-300 MHz. In step 404, a first distancebetween the source and detector is chosen. In step 406, measurements ofDC and AC intensity and phase shift are made at the first distance andfrequency. In step 408, the distance is altered and measurements arerepeated at the same modulation frequency. In step 408A, the criteriaparameter, Q, is calculated. If Q has a value near 1, then scatteringand absorption coefficients are calculated from the measurements and alinear analysis of the measurements at different distances. Measurementsare then repeated at a different modulation frequency in step 410. Instep 412, calculations are performed for the different modulationfrequencies. An average value of absorption coefficient μ_(s)′ at thewavelength of light used is then calculated. Measurements may berepeated using light of a different wavelengths.

Equations (10a)-(10d) also show that at fixed modulation frequency, ω,the values of${\ln\left( {\frac{r}{r_{0}}\quad{DC}_{rel}} \right)},\quad{\ln\left( {\frac{r}{r_{0}}\quad{AC}_{rel}} \right)},$ln(Mod,_(rel))and (PS_(rel)) are linear functions of the difference ofdistance between detectors, (r−r₀). From a plot of these values,${\ln\left( {\frac{r}{r_{0}}\quad{DC}_{rel}} \right)},\quad{\ln\left( {\frac{r}{r_{0}}\quad{AC}_{rel}} \right)},$ln(Mod_(rel))and (PS_(rel)), as a function of (r−r₀), the slopes can bedetermined from linear regression. We term these slopes, K_(DC), K_(AC),K_(Mod) and K_(PS), respectively. These slopes are predictive of theoptical properties, using the following equations: $\begin{matrix}{K_{D\quad C} = {- \sqrt{3{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}}}} & \left( {10a} \right) \\{K_{A\quad C} = {- {\sqrt{\frac{3}{2}{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}}\left\lbrack {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} + 1} \right\rbrack}^{1/2}}} & \left( {10b} \right) \\{K_{Mod} = {- {\sqrt{\frac{3}{2}{\mu_{a}\left( {\mu_{a} + \mu_{a}^{\prime}} \right)}}\quad\left\lbrack {\sqrt{2} - \left( {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} + 1} \right)^{1/2}} \right\rbrack}}} & \left( {10c} \right) \\{K_{PS} = {\sqrt{\frac{3}{2}{\mu_{a}\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}}\left\lbrack {\sqrt{1 + \left( \frac{\omega}{v\quad\mu_{a}} \right)^{2}} - 1} \right\rbrack}^{1/2}} & \left( {10d} \right)\end{matrix}$

If the refraction index of the medium is known, the absorptioncoefficient, μ_(a), and isotropic scattering coefficient, μ_(s)′, can bedetermined given the modulation frequency, ω. The advantage of the MDmethod is that an analytical solution of optical parameters can bedirectly derived from the above equations and nonlinear regression canbe avoided.

The above four equations contain only two unknowns, and therefore, theyare not independent. The following equation can be obtained fromEquation (10) for checking the consistency of FDPM measurements:$\begin{matrix}{\frac{K_{AC}^{2} - K_{PS}^{2}}{K_{DC}^{2}} = {K = 1}} & (11)\end{matrix}$We term the criteria parameter for the MD method, K.

In practice, the values of μ_(a) and μ_(s)′ can be obtained using fourdifferent linear regression strategies: simultaneous regression of (1)DC+AC; (2) DC+PS; (3) AC+PS; and (4) Mod+PS experimental data. Fromappropriate algebra transformation, the analytical solutions of opticalparameters can be expressed as: $\begin{matrix}{{AC} - {{DC}:}} & \quad \\{\mu_{a} = {{\frac{\omega\quad n}{c}\left\{ {\left\lbrack {{2\left( {K_{AC}/K_{DC}} \right)^{2}} - 1} \right\rbrack^{2} - 1} \right\}^{1/2}\quad\mu_{s}^{\prime}} = {\frac{K_{DC}^{2}}{3\mu_{a}} - \mu_{a}}}} & \left( {12a} \right) \\{{DC} - {{PS}:}} & \quad \\{\mu_{a} = {{\frac{\omega\quad n}{c}\left\{ {\left\lbrack {{2\left( {K_{PS}/K_{DC}} \right)^{2}} + 1} \right\rbrack^{2} - 1} \right\}^{1/2}\quad\mu_{s}^{\prime}} = {\frac{K_{DC}^{2}}{3\mu_{a}} - \mu_{a}}}} & \left( {12b} \right) \\{{AC} - {{PS}:}} & \quad \\{\mu_{a} = {{{\frac{\omega\quad n}{c}\left\lbrack {\left( \frac{K_{AC}^{2} + K_{PS}^{2}}{K_{AC}^{2} - K_{PS}^{2}} \right)^{2} - 1} \right\rbrack}^{{- 1}/2}\quad\mu_{s}^{\prime}} = {\frac{K_{AC}^{2} - K_{PS}^{2}}{3\mu_{a}} - \mu_{a}}}} & \left( {12c} \right)\end{matrix}$Where n is the refraction index of the scattering medium and c is thelight speed in a vacuum.

For Mod and PS data, the analytical solution for the optical propertiescannot be directly obtained owing to the lack of linear dependence onμ_(a): $\begin{matrix}{\frac{K_{Mod}}{K_{PS}} = {{\frac{\sqrt{2} - \sqrt{T + 1}}{\sqrt{T - 1}}\quad{where}\quad T} = \sqrt{1 + \left( \frac{\omega\quad n}{c\quad\mu_{a}} \right)^{2}}}} & \left( {12d} \right)\end{matrix}$However, the numerical solution of μ₁ can be easily obtained, and μ_(s)′can be subsequently obtained from equation (10d). Values of μ_(a) andμ_(s)′ can also be obtained from simultaneously regressing all AC, DC,and PS data by the least squares method.

FIG. 5 illustrates use of the criteria parameter, Q, as a function ofrelative distance r−r₀ and modulation frequency ω for a typical FDPMmeasurement. Three experimental points at small relative distances r−r₀and small frequencies ω are out of the range of this figure because thevalues of criteria parameter Q at these points are large. In thisfigure, the values of function Q are close to one for most of theexperimental data with the exception of those data at small relativedistances that are outside preferred operating conditions, and whichsatisfy the assumptions of the photon diffusion model. FIG. 5illustrates typical experimental results. For each experiment, analysisof the criteria parameter provides information of: (1) whether abnormalmeasurement error exists during FDPM experiment; (2) which ranges ofmodulation frequency and relative distance are suitable for FDPMexperiments with the particular sample; and (3) what parts of themeasurement data can be used to give accurate and reliable opticalproperties. So Eq, 9 provides an important criterion to check theinstrument performance, to avoid system and abnormal measurement error,and to find optimal operation conditions during FDPM measurements.

Compared with theoretical predictions, simultaneous regression of DC+PS,AC+PS or DC+AC+PS data can provide good derived optical parameters athigh modulation frequency and large relative distances for both the MFand MD methods. In order to increase the signal-to-noise ratio, and toimprove the accuracy and precision of parameter estimation, the MD andMF methods may be combined. For a MF measurement, the experiment can beperformed at several different relative distances (called the combinedMF method); while for a MD measurement, the experiment can be performedat several different modulation frequencies (called the combined MDmethod). The averaged accuracy and precision of estimated scatteringcoefficients using combined MF and MD measurements with differentregression approaches and relative distances ranging from 2.5 mm to 8 mmfor MF measurements and modulation frequencies ranging from 60 MHz to100 MHz for MD measurements showed that these ranges avoid large noisesarising at low frequencies and small relative distances (“PreciseAnalysis of Frequency Domain Photon Migration Measurements forCharacterization of Concentrated Colloidal Suspensions,” Rev. Sci.Instrum., 73(2): 383-393, 2002). Compared with MD and MF methods usedalone, all these combined approaches can improve the accuracy andprecision of estimated isotropic scattering coefficients. The averagedaccuracy of the combined MD and MF methods are comparable to theuncertainty of the calculated μ_(s)′ values by Mie theory (1.44%).Compared to the combined MF method, the combined MD method hasadvantages that it avoids the non-linear regression. Furthermore, linearplot and criteria parameter K (Eq. 11) can easily be used to checkexperiment conditions and consistency of experimental data with the MDmethod. Since simultaneous regression of DC+PS data or AC+PS data ofcombined MD measurements provides estimation of uncertainties ofisotropic scattering coefficient by standard error analysis of a linearmodel, both regression strategies are preferred for data analysis.

When the final value of scattering coefficient of a sample is determinedat a selected wavelength of the incident light, that value is insertedin equations that model the scattering properties of the colloidalsample to be investigated. The scattering properties depend on volumeconcentration of a disperse phase, the particle size distribution of theparticles and, at higher volume concentrations, the interactions betweenthe particles.

In past methods, the particle size distribution (PSD) of a sample hasbeen assumed a priori to be either normal, log-normal or Weibull. Theparameters that describe each of these distributions were determined bya best-fit procedure using non-linear least-squares regression. Usingthe following procedure, the particle size distribution of a sample canbe determined without an a priori assumption.

For non-interacting suspension (i.e., for S=1), upon measuring isotropicscattering coefficients at different wavelengths using FDPM, the inverseproblem of determining f(x) and φ requires solution of the integralequation (4). Eq. (4) can be generalized as the following set ofintegral equations,μ_(s)′(λ_(i))=∫_(a) ₁ ^(a) ² K(x,λ _(i))p(x)dx+e _(i) i=1,2, . . . M  (13)where K(x, λ_(i) is the kernel function, which is given by:$\begin{matrix}{{{K\left( {x,\lambda_{i}} \right)} = {{{\frac{3{Q_{scat}\left( {\lambda_{i},x} \right)}}{2x}\left\lbrack {1 - {g\left( {\lambda_{i},x} \right)}} \right\rbrack}\quad i} = 1}},\ldots\quad,M} & (14)\end{matrix}$Here μ_(s)′ (λ_(i) is the discrete set of measurement data at Mdifferent wavelengths; p(x) is an unknown function to be determined; a₁and a₂ are the size limits within which the size distribution lies; ande_(i) is the measurement error.

Eqs (13) and (14) constitute the inversion problem for recovering anunknown function p(x) from discrete measurement data μ_(s)′(λ_(i)). Inthe method disclosed herein, the unknown function p(x) is the product,φf(x), instead of f(x) and φ individually, since f(x) and φ can beseparated by the normalized distribution constraint:∫φf(x)dx=φ  (15)The unknown function p(x) is approximated by a linear combination ofB-spline base functions, B_(j) ^(m)(x), with coefficients, c_(j), whichare given by: $\begin{matrix}{{p(x)} = {{\sum\limits_{j = 1}^{N}{c_{j}{B_{j}^{m}(x)}}} + {ɛ(x)}}} & (16)\end{matrix}$Here ε(x) is the approximation error, and m, N are the order anddimension of the B-spline respectively. Since B-spline can accuratelyapproximate any smooth function, the use of B-spline functions reducesthe number of independent variables to be determined and stabilizes theinversion process.

Using the B-spline approximation, Eq. 13 can be written as:$\begin{matrix}{{{\mu_{s}^{\prime}\left( \lambda_{i} \right)} = {{{\sum\limits_{j = 1}^{N}{c_{j}{\int{{K_{i}\left( {x,\lambda_{i}} \right)}{B_{i}^{m}(x)}{\mathbb{d}x}}}}} + {E_{i}\quad i}} = 1}},\ldots\quad,M} & (17)\end{matrix}$where E_(i) is the summation of measurement-error and model(approximation) error. In matrix notation, equation (15) is given by:b=Af+E   (18)where b is the measurement data and b_(i)=μ_(s)′(λ_(i)); f is an unknownfunction and f_(j)=c_(j); E is error; and A is the response model, whosecomponents are:a _(i,j) =∫K _(i)(x,λ _(i))B _(j) ^(m)(x)dx   (19)In general, the solution to equation (18) can be formulated as aclassical least squares problem:min{∥Af−b∥₂ ²}  (20)However, since direct solution of equation (19) can yield an unstablesolution when measurement data is contaminated with noise,regularization is necessary for reducing the influence of the noise inthe inversion process.

To stabilize the inverse solution, we employ Tikhonov regularization inwhich the regularized solution fα is defined as the solution of thefollowing optimizing problem:min{∥Af−b∥₂ ²+α²∥Lf∥₂ ²}  (21)Here L is an appropriately chosen regularization matrix. Theregularization parameter a controls the weight given to minimization ofthe regularization term, ∥Lf_(reg)∥₂ relative to minimization of theresidual term, ∥Af_(reg)−b∥₂.

Different ways to choose the regularization matrix L result in differentregularization methods, such as zeroth, first, and second order, andentropy regularization. In the method disclosed herein, the second orderregularization form is chosen since it provides the smoothness of thesolution, which is often called a smoothing constraint: $\begin{matrix}{{{Lf}_{reg}}_{2}^{2} = {\sum\limits_{j = 3}^{N}\left( {f_{j} - {2f_{j - 1}} + f_{j - 2}} \right)^{2}}} & (22)\end{matrix}$However, the smoothing constraint alone cannot guarantee stable,non-negative solutions of the size distribution. To avoid a negativesolution, the non-negativity constraint should be included:$\begin{matrix}{{p(x)} = {{\sum\limits_{j = 1}^{N}{f_{j}{B_{j}^{m}(x)}}} \geq 0}} & (23)\end{matrix}$

Equations (21), (22) and (23) constitute the inverse problem in whichcoefficient f of the B-spline expansion needs to be determined. Fromequations (15) and (16), the PSD and volume fraction are subsequentlyobtained. If the regularization parameter α is known, the inverseproblem is a linear least squares problem with linear inequality. Theactive-set quadratic programming (QP) method may be used to solve thePSD inverse problem. The success of this regularization method dependson the choice of optimal regularization parameter α.

The choice of optimal regularization parameter is important to solve theill-posed inversion problem. On one hand, if the parameter is too small,the noise cannot be filtered out and spurious oscillations in thesolution exist. On the other hand, if the parameter is too large, somedata are filtered out as well as noise and the solution isover-smoothed.

If information about error norm ∥E∥₂ (including model and experimentalerrors) is available, the discrepancy principle can be used to chooseregularization parameters, such that the residual norm for theregularized solution equates the error norm. However, this techniqueclearly depends on the ability to correctly estimate errors in themeasurement data. Generalized cross validation (GCV) is a popular methodfor choosing the regularization parameter since it does not require anestimation of the error norm. The GCV method is based on statisticalconsiderations that a good regularization parameter can predict missingdata values. However, GCV cannot be expected to perform well in the casewhere there are relatively few measurements and the GCV function mayhave a flat minimum. Consequently it may be difficult to locate thevalue of a global minimum.

The L-curve method is another parameter-choosing method that isindependent of error estimates. The L-curve is a plot of regularizationterm ∥Lf_(reg)∥₂ versus residual term ∥Af_(reg)−b∥₂ for allregularization parameters. As shown in FIG. 6, an optimal solution liesnear the corner of the “L-curve” since the solution strikes a balancebetween the solution smoothness and the residual error in this region.In most applications, the L-curve method is used for regularizationmethods which employ a smoothing constraint only. In the methoddisclosed herein, the L-curve method is used to choose theregularization parameter for the PSD inversion problem contained by bothsmoothness and non-negativity constraints. Results of numericalsimulation to investigate the performance of the method are reported inthe paper “Inversion Algorithms for Particle Sizing with PhotonMigration Measurement,” AIChE J., July, 2001, Vol. 47, No. 7, pp.1487-1498, which is incorporated by reference herein.

FIG. 7 shows a flowchart of the method for determining volume fractionand particle size distribution for a sample of non-interacting particleshaving any size distribution—single modal or multi-modal. If it issuspected that the particles in a sample may not have a simple particlesize distribution, i.e., may be polydisperse or multi-modal, theprocedure of FIG. 7 should be used in analysis of the scattering data.The method is described in the paper “Inversion Algorithms for ParticleSizing with Photon Migration Measurement,” AIChE J. Vol. 47, No. 7, 1487(July 2001).

Referring to FIG. 7, in step 502, values of refractive indices,wavelength of the incident light, the measured value of scatteringcoefficient and an estimate of the upper and lower bounds of particlessizes in the sample are entered. In step 504, an approximation of theunknown function p(x) is expressed as a linear combination of B-splinebase functions with coefficients c. The value of N is to be no greaterthan the number of wavelength measurements. In step 506, kernel functionK is calculated using Mie theory, which is well known in the art.Matrices A and b are then formulated as shown in the figure. In step 508an optimization procedure employing a regularization matrix L is usedwith second order regularization. Then in step 510, the optimizationproblem is solved by an active-set quadratic programming method, asshown. In step 512, an L-curve or another approach is used to choose thebest regularization parameter and to determine the solution vector C. Instep 514 the value of p(x) is determined, and in step 516 the value ofvolume fraction is determined by integrating the function p(x).

Simulation studies showed that: (1) the B-spline expansions provide agood approximation of size distributions; (2) the inversion algorithmcan accurately recover the PSD and volume fraction if there is nomeasurement error; and (3) the effects of model-errors can be neglected,reducing the need for regularization significantly. The use of theB-spline on actual experimental data is also demonstrated in the paper“Inversion Algorithms for Particle Sizing with Photon MigrationMeasurement,” AIChE J., July, 2001, Vol. 47, No. 7, pp.1487-1498.

The inverse algorithm can effectively suppress noise propagation, andprovide reasonable inversion results for normal and log-normaldistributions within as high as 10% measurement error in isotropicscattering coefficients. As the noise level increases, the sizedistributions are over-smoothed, which may be due to the L-curve method.Although the accuracy of recovered size distributions decreases as thenoise level increases, the recovered volume fractions are accurate, andthe moments of recovered PSD (i.e., mean size and standard deviation)are close to the true values.

From solutions using synthetic measurements, both the smooth andnon-negativity constraints are important for PSD inversion from FDPMmeasurements. The original Tikhonov regularization method does notinclude the non-negativity constraints. When used in the inversesolution for FDPM measurements, it produces unrealistic negativesolutions, especially at the limits of particle size distribution.

As described in the paper “Inversion Algorithms for Particle Sizing withPhoton Migration Measurement,” AIChE J., July, 2001, Vol. 47, No. 7, pp.1487-1498, two Gaussian size-distributions with the same mean sizes butdifferent standard deviations, using synthetic measurements, showed thatat noise levels of 1%, 5% and 10%, the effects of polydispersity on thequality of inversion are small, and the method described herein canrecover both narrow and broad distributions at different noise levels.Also, the method disclosed herein can recover a bimodalsize-distribution. The recovered results were in good agreement withassumed true values for the bimodal size distribution and the modelerrors were negligible. Furthermore the results demonstrate that the apriori assumption of a PSD can be avoided.

A study of the effects of measurement errors for bimodal sizedistribution showed that at low noise level (less than 5%), the firstmode shows good agreement with the true value, while the second mode isfar away from the true value. Compared to the inversion of unimodalnormal distribution, the effects of noise are larger for recovery ofbimodal normal distributions, which means that recovery of multi-modaldistribution is more sensitive to measurement error than recovery of aunimodal distribution. Hence, more accurate and more numerousmeasurement data are required in order to recover a multi-modal sizedistribution. Although agreement between retrieved PSD and assumed PSDdecreases as noise level increases, the two modes of distribution can besuccessfully reconstructed with random measurement error as high as 10%.

Examples of recovery of PSDs for both synthetic data and experimentaldata are provided in the paper “Inversion Algorithms for Particle Sizingwith Photon Migration Measurement,” AIChE J., July, 2001, Vol. 47, No.7, pp. 1487-1498. Particle-sizing results for experimental data ontitanium dioxide samples are compared with results of previous methodsusing dynamic light scattering and X-ray measurements.

The procedure described in FIG. 7 is applicable to samples at volumeconcentrations low enough such that particles are non-interacting. Whenparticles are interacting, the nature of interaction forces between theparticles can also be investigated by FDPM measurements using Eq. 5, butthe partial structure factor, S_(i,j), included in the equation, whichis a direct measure of the local ordering of colloidal particles, willbecome less than unity.

For a polydisperse system, since particle interactions occur betweencolloids not only with the same particle size but also with differentparticle sizes, the partial structure factor, S_(i,j), should beconsidered, where i and j represent particles with different sizes.

Although the form factor, F_(i,j), in Eq. 5 can be calculated byclassical Mie scattering theory, it is difficult to extract the partialstructure factor, S_(i,j), from Eq. (5) using the measured isotropicscattering coefficient, μ_(s)′, from FDPM measurements at severaldifferent wavelengths, owing to ill-posedness of the Fredholm integralequation. Instead, the angle-integrated structure factor, <S(q)>, can bedefined from the ratio of scattering coefficients computed or measuredwhen particle interactions are accounted for (i.e.,μ_(s,interacting)′)and when they are not (i.e.,μ_(s,non-int eracting)′): $\begin{matrix}{{< {S(q)}>=\frac{\left( \mu_{s}^{\prime} \right)_{interacting}}{\left( \mu_{s}^{\prime} \right)_{{non} - {interacting}}}} = \frac{\begin{matrix}{\int_{0}^{\pi}{\sin\quad{\theta\left( {1 - {\cos\quad\theta}} \right)}{\int_{0}^{\infty}{\frac{f\left( x_{i} \right)}{x_{i}^{3}}{\int_{0}^{\infty}\frac{f\left( x_{j} \right)}{x_{j}^{3}}}}}}} \\{{F_{i,j}\left( {n,x_{i},x_{j},\lambda,\theta} \right)} \cdot} \\{{S_{i,j}\left( {n,x_{i},x_{j},\lambda,\theta} \right)}{\mathbb{d}x_{j}}{\mathbb{d}x_{i}}{\mathbb{d}\theta}}\end{matrix}}{\begin{matrix}{\int_{0}^{\infty}{\frac{f(x)}{x^{3}}{\int_{0}^{\pi}{F\left( {n,x,\lambda,\theta} \right)}}}} \\{\sin\quad{\theta\left( {1 - {\cos\quad\theta}} \right)}{\mathbb{d}\quad\theta}{\mathbb{d}x}}\end{matrix}}} & (24)\end{matrix}$Hence, the angle-integrated structure factor <S(q)> can beexperimentally obtained using FDPM-measured μ_(s)′ of concentratedsuspensions as well as from theoretical prediction from a model ofparticle interaction or partial structure factor, S_(i,j). By comparingangle-integrated structure factors obtained experimentally with thatthose theoretically predicted, one can examine the model of S_(i,j) toinvestigate particle-particle interactions in colloidal suspensions. Wefound that particle interactions arise primarily due top the excludedvolume interaction at high volume fractions, and the polydisperse hardsphere Percus-Yevick (HSPY) model is suitable for accounting forparticle interactions in dense, polydisperse suspensions.

For an interacting system, if the form of PSD is known (e.g., Gaussian,log normal, etc.), we can recover the PSD and volume fraction byregressing the unknown parameters of PSD using nonlinear optimizationmethods. The algorithm to solve this inverse problem involves theminimization of a merit function, χ²: $\begin{matrix}{\chi^{2} = {\sum\limits_{j = 1}^{M}\left\lbrack {\left( \mu_{s}^{\prime} \right)_{j}^{O} - \left( \mu_{s}^{\prime} \right)_{j}^{C}} \right\rbrack^{2}}} & (25)\end{matrix}$where (μ_(s)′)_(j) ^(o) is the experimentally observed scatteringcoefficient and (μ_(s)′)_(j) ^(c): is the value calculated from equation(5), based on the current guess for size distribution f(x).

If the form of PSD is unknown, for interacting systems, Eq. (5) is to besolved for recovery of volume fraction φ and particle size distribution,f(x). Since S_(i,j) is a function of unknown φ, Equation (5) is anonlinear Fredholm integral equation. A general method for solving thenon-linear inverse problem for any arbitrary model is illustrated inFIG. 8. In step 701, data and the measured value of the scatteringcoefficient at a selected wavelength are input. In step 702, values of φand f(x) are assumed. In step 703, the nonlinear inverse problem iscalculated to recover values of f(x) and φ′. In step 704, recoveredvalues are compared with assumed values. If the difference is not lessthan a selected value, update values are then selected in step 705. Ifthe difference is less than the selected value, values of f(x) and φ areoutput in step 706.

In order to avoid solving a nonlinear inverse problem directly, severalapproximations, such as the Decoupling Approximation (DA) and the LocalMonodisperse Approximation (LMA), can be used to simplify thecomplicated polydisperse Hard Sphere Percus-Yevick (HSPY) model toaccount for volume exclusion. In the DA model, it is assumed that thepositions of the particles are independent of their size so that thepartial structure factor is given by:S _(i,j)(n,λ,x _(i) ,x _(j),φ)=S(n,λ,{overscore (x)},φ)   (26)where {overscore (x)} is the average diameter of all colloids in thepolydisperse system. For the LMA model, it is assumed that the positionsof particles are completely correlated with their sizes so that eachparticle is surrounded by particles with the same size. The partialstructure factor in the LMA model can therefore be written as:$\begin{matrix}{{S_{i,j}\left( {n,\lambda,x_{i},x_{j},\phi} \right)} = \left\{ \begin{matrix}{S_{i}\left( {n,\lambda,x_{i},\phi} \right)} & {i = j} \\0 & {i \neq j}\end{matrix} \right.} & (27)\end{matrix}$In both DA and LMA models, the effective structure factor can becalculated by the monodisperse HSPY model. Correspondingly, Equation (5)can be written as: $\begin{matrix}\begin{matrix}{{\mu_{s}^{\prime}(\lambda)} = {\int_{0}^{\infty}{\frac{12\quad\phi\quad{f(x)}}{k^{2}x^{3}}\left\lbrack {\int_{0}^{\pi}{{F\left( {n,x,\lambda,\theta} \right)}{S\left( {n,x_{eff},\lambda,\phi_{eff}} \right)}\sin\quad\theta}} \right.}}} \\{\left. {\left( {1 - {\cos\quad\theta}} \right){\mathbb{d}\theta}} \right\rbrack{\mathbb{d}x}}\end{matrix} & (28)\end{matrix}$where x_(eff) and φ_(eff) are the effective diameter and volumefractions that arise due to particle interactions, andφ=φ_(eff)/(x_(eff)/x)³. For the DA model, we assume thatx_(eff)={overscore (x)}, since the approximation simply states that allinteractions are approximated by the interactions of particles of themean size of the distribution. For the LMA model, we assume: x_(eff)=x,and that only like particles of the same size interact. Hence the volumefraction φ_(eff) is essentially the volume fraction of particles withx_(eff).

The DA and LMA models have been examined by FDPM measurements in densesurfactant-stabilized suspensions of polystyrene, and examples of modelcomparisons are shown in FIG. 11, from the paper “Approach for ParticleSizing in Dense Polydisperse Colloidal Suspension Using MultipleScattered Light,” Langmuir 2001, 17, 6142-6147, which is incorporated byreference herein. The figure shows that, at volume fractions <0.1 the DAmodel is a good approximation to the data and the HSPY model. Similarresults can be inferred for the LMA model if the volume fraction is lessthan 0.2. The mismatch between predictions by the DA and the LMA modelsand experimental data at high volume fractions is due to the fact thatthe effective diameter of the particles, because of particleinteractions, cannot be assumed to equal their physical diameter.Correspondingly, if we introduce the interaction parameter, which isdefined as β=x_(eff)/{overscore (x)} for the DA model, or β=x_(eff)/xfor the LMA model, the accuracy of the DA and the LMA models can besignificantly improved. Table 1 illustrates the interaction parameters,β, of the DA and LMA models by fitting FDPM experimental data to themodels with known volume fraction and PSD, where experimental datameasured at wavelengths of 687 nm and 828 nm were obtained from the samesamples. We can see that, at each volume fraction, the interactionparameters fitted from independent data acquired at two differentwavelengths are in close agreement. This indicates that, by introducingthe interaction parameter, the DA and LMA model can accurately predictthe FDPM experimental data even at high volume fractions (>0.25).

A flowchart for determining φ, f(x) and β using the LMA model is shownin FIG. 9(a) and FIG. 9(b). The parameter β is set equal to 1 and instep 902 an initial value of φ is guessed and set equal to φ_(eff). Instep 904, input physical data, measured data and assumption are input asshown. In step 906, a B-spline approximation of p(x) of dimension N isintroduced. In step 908, the kernel function from Mie theory iscalculated and S is calculated from the LMA model, using matrices asshown in the figure. In step 910, the optimization approach is used withuse of matrix L and second order regularization. In step 912 theoptimization problem is solved using the equations shown and thenon-negative constraint of p(x). In step 914 the L-curve or an alternateapproach is used to choose the best regularization parameter. In step916, p(x) is calculated. In step 918, the value of resulting φ iscalculated, called φ′. In step 920, the calculated value of φ iscompared with the initial assumed value to determine if the differenceis less than the convergence error. If it is not, in step 920A a newvalue of φ is chosen and a new value of φ_(eff) is calculated in step920B. If the difference is less than convergence error, the value of β′is calculated in step 922. The check for convergence of β is determinedin step 924. If the difference is greater than convergence error, a newvalue of β is chosen in step 926. If it is not, then the final values off(x), φ and β have been determined.

Similarly, a flowchart for determining φ, f(x) and β using the DA modelis shown in FIG. 10. The parameter β is set equal to 1 and in step 1002an initial value of φ and a mean size (x bar) are guessed. In step 1004,input physical data, measured data and assumptions are input as shown.In step 1006, a B-spline approximation of p(x) of dimension N isintroduced. In step 1008, the kernel function F from Mie theory iscalculated and S is calculated from the DA model, using matrices asshown in the is figure. In step 1010, the optimization approach is usedwith use of matrix L and second order regularization. In step 1012 theoptimization problem is solved using the equations shown and thenon-negative constraint of p(x). In step 1014 the L-curve or analternate approach is used to chose the best regularization parameterand the solution vector. In step 1016, p(x) is calculated. In step 1018,the value of resulting φ is calculated, called φ′. Also, (x bar) iscalculated from f(x). In step 1020, the calculated value of φ iscompared with the initial assumed value to determine is the differenceis less than the convergence error. If it is not, in step 1020A a newvalue of φ and/or (x bar)is chosen and input at step 1006. If thedifference is less than convergence error, the value of β′ is calculatedin step 1022. The check for convergence of β is determined in step 1024.If the difference is greater than convergence error, a new value of β ischosen in step 1026. If it is not, then the final values of f(x), φ andβ have been determined. TABLE 1 Interaction Interaction parameterparameter Volume (DA model) (LMA model) fraction 687 nm 828 nm 687 nm828 nm 0.0510 1.127 1.069 1.161 1.096 0.1007 1.017 1.014 1.048 1.0380.1516 0.994 0.992 1.021 1.013 0.2017 0.976 0.978 1.002 0.996 0.25360.964 0.963 0.988 0.979 0.3015 0.958 0.958 0.979 0.972 0.3531 0.9530.954 0.971 0.966 0.4032 0.949 0.949 0.965 0.959

EXAMPLE

A polydisperse sample of concentrated, surfactant stabilized,polystyrene suspension was provided by Dow Chemical, Midland, Mich., andused as received. The mean particle size, deviation and polydispersityof the polystyrene sample were measured by using DLS (ZETASIZER 3000,Malvern Instruments, U.K.) at dilute suspensions (˜0.01%), and theresults are shown in Table 2. TABLE 2 Mean size Deviation (nm) (nm)Polydispersity (%) DLS results 143.5 ± 0.5 24.0 ± 1.2 16.7 ± 0.8

FDPM measurements were conducted at various particle concentrationsranging from 5% to 40% solids by volume. The stock solutions werediluted by deionized ultra filtered (DUF) water. The exact volumefractions of samples were determined by evaporation measurements, andare shown in Table 3. TABLE 3 Crude volume   5%   10%   15%   20%   25%  30% fraction Exact volume 5.10% 10.07% 15.16% 20.17% 25.36% 30.15%fraction

Based on the photon diffusion model, several different experimentalapproaches and regression strategies were used to determine theabsorption coefficient, μ_(a), and isotropic scattering coefficient,μ_(x)′. In the experiments, multiple distance (MD) measurements wereperformed at twelve different source-detector distances, ranging from 7mm to 11 mm. In order to improve the precision and accuracy, MDmeasurements were performed at several different modulation frequenciesranging from 50 to 100 MHz. The average of FDPM-determined opticalparameters and their uncertainties were calculated in order to decreasemeasurement errors. This method may be called the combined MDmeasurement, and gives the highest accuracy and precision of isotropicscattering coefficients.

FIG. 12 illustrates a plot of isotropic scattering coefficient measuredby FDPM measurements as a function of volume fraction φ for the latexsamples at two different wavelengths.

Besides the isotropic scattering coefficients, their uncertainties canalso be derived from FDPM experimental data. FIG. 13 illustrates theprecision of the FDPM-measured scattering coefficients at twowavelengths and a series of volume fractions. Here, the precision isdefined as relative standard deviation of regressed isotropic scatteringcoefficients. The results show that the precision of FDPM-measuredisotropic scattering coefficients is less than 0.4% at both wavelengthsin all the measurements except at the lowest volume fraction.

Using FDPM-measured isotropic scattering coefficients, theangle-integrated structure factor, <S(q)>, was experimentally determinedusing Eq. 24. FIG. 14 illustrates a plot of <S(q)> as a function ofvolume fraction at two wavelengths for latex samples. The theoreticalpredictions using the polydisperse HSPY model, also shown in FIG. 14,indicate that the HSPY model represents a good approximation fordescription of particle interactions. Although small differences areobserved between the experimental results and theoretical predictions,our results show that: (1) angle-integrated structure factors decreasewith increasing volume fractions at the same wavelength due to the factthat particle interactions become more significant with increasingvolume fraction; and (2) angle-integrated structure factors decreasewith increasing wavelength at the same volume fraction due to the factthat as wavelength decreases, wave vector increases and <S(q)> becomelarger. The small differences between the experimental results andtheoretical predictions may be due to: (1) the existence of other typesof particle interactions; and (2) possible discrepancies of DLS-measuredmean size of latex particles used for calculation of μ_(s)′, which wouldnot be accounted for in this analysis. Nonetheless, FIG. 14 shows thatparticle interactions arise primarily due to the excluded volumeinteraction at high volume fractions, and the polydisperse HSPY model issuitable for accounting for particle interactions in this polydisperselatex sample.

Since the polydisperse HSPY model is a good approximation for modelingparticle interactions of colloidal suspensions at high volume fractions,we can incorporate this model into the expression of isotropicscattering coefficient (i.e., Eq. 5) to recover an unknown particle sizedistribution. In this work, we recovered the particle size distributionof latex samples by fitting FDPM-measured μ_(s)′ at various volumefractions at a wavelength of 687 nm or 828 nm using nonlinearregression. FIG. 15 illustrates the FDPM-measured and fitted isotropicscattering coefficients using the polydisperse HSPY model as a functionof volume fraction at wavelengths of 687 nm. For comparison, the dottedline represents the theoretical value of the isotropic scatteringcoefficient neglecting particle interactions (i.e., S(q)=1), aspredicted by using Eq. 5. It can be seen that particle interactionssignificantly contribute to the FDPM-measured isotropic scatteringcoefficient at high volume fraction, leading to a large deviation of themeasured scattering coefficient from prediction by Mie theory. Table 4lists the recovered mean sizes and deviations from FDPM measurements atwavelengths of 687 nm and 828 nm, which agree well with the DLS resultsmeasured at dilute latex sample (˜0.01% volume fraction). Theseinversion results demonstrate that the FDPM method with the polydisperseHSPY structure factor model accounting for dominant particleinteractions can be successfully used for recovery of particle sizedistribution of interacting colloidal suspensions at high volumefractions. TABLE 4 Mean size Deviation Polydispersity (nm) (nm) (%) FDPMresults (687 nm) 142.23 25.46 17.9 FDPM results (828 nm) 143.14 21.8415.3

The ability to recover particle size distribution of colloidalsuspensions at high volume fractions depends largely upon finding anappropriate model for considering effects of particle-particleinteractions. Most of available particle sizing methods requireempirical calibration when sizing at high volume fractions due tounavailable models for considering particle interactions. Here it isdemonstrated that the FDPM technique is a useful tool for probingparticle interactions in polydisperse, concentrated, multiple scatteringcolloidal suspensions. Compared with theoretical predictions, thepolydisperse HSPY model is suitable for description of particleinteractions and local microstructure of colloidal suspensions at highvolume fractions. By incorporating the polydisperse HSPY model toaccount for particle interactions, particle size distribution from FDPMmeasurements at high volume fractions (5-40%) were determined. Theresults agree well with DLS measurement at low volume fractions(˜0.01%). While the previous work demonstrated that the FDPM techniquecan be used for particle sizing of multiple-scattering yetnon-interacting colloidal suspension (<5% volume fraction), this workextends FDPM techniques to characterize interacting suspensions at highvolume fraction. Compared with other methods, the FDPM technique is mostsuitable for particle sizing of colloidal suspensions at high volumefractions, due to the fact that: (1) it is based upon first principles;(2) it considers the effects of both multiple scattering and particleinteractions; and (3) it is self-calibrating, and no externalcalibration is required.

In the above discussion, the particle interactions that were consideredwere brought about by volume exclusion effects; that is, by the factthat no two particles can occupy the same volume. As a result of thisrestriction, the particles in suspension become “ordered” or structuredand the structure factor, S, or the partial structure factors, S_(ij),become less than one. For most systems, the primary order in asuspension is caused by the volume exclusion effect. However, other moresubtle interactions which arise due to interparticle forces can furtheradd to the “order” and structure. The first-principles models forcomputing S and S_(ij) for neutron and x-ray scattering are available orare becoming available in the open literature. Interaction forces thatarise owing to electrostatics, depletion forces, van der Waals forces,etc., impact the order and, therefore, the structure.

FDPM offers an opportunity to assess these interactions and theparameters that govern them. When a particle size is known or determinedby an auxiliary method, the scattering within a dense suspension can beaccurately measured with FDPM and used to assess the forces which existin dense suspensions. In the following, we demonstrate the approach byprobing the effective surface charge of a dense suspension of negativelycharged polystyrene.

FIG. 16 illustrates the effect of salt concentration of a polystyrenesample on the isotropic scattering coefficient. At high volumefractions, scattering in the dialyzed sample is decreased because ofgreater interaction between the particles, which is caused byelectrostatic repulsion and volume exclusion. When the saltconcentration increases, this “screens” the electrostatic interactionand results in a reduction of order and consequently an increase in thescattering coefficient. The following discussion includes thecharacterization of dense suspensions under repulsive electrostaticinteractions between particles. The discussion also applies toapplication of FDPM to assess changes in scattering due to electrostaticinteractions, depletion forces, van der Waals force, steric forces andinteraction parameters that relate to these forces. The interactionsamong colloidal particles decide the local structure, which can impactrheology, colloidal stability, and therefore the final performance ofthe industrial colloidal products. When characterizing colloidal systemsusing scattering techniques, such as x-ray and light scattering, thestructure of the suspension has strong influence over scatteringefficiency, especially at condensed cases.

The particle-particle interaction of charged particles can be consideredas the addition of two parts: hard sphere interaction and electrostaticor other force interaction. Electrostatic interaction comes from theoverlap of electrical double layers of particles while hard sphereinteraction comes from volume exclusion effects or the simple fact thattwo particles cannot occupy the same volume. Upon evaluating a firstprinciples model to account for the force, a structure factor or partialstructure factor can be determined which is a function of the parametergoverning the interparticle forces. In this example, the parameterimpacting electrostatic interaction is the effective surface charge. Thefollowing outlines how the structure factor can be used to reflect theseattractive and repulsive interactions from first principles.

Consider two unlike, interacting particles, i and j of diameter x_(i)and x_(j). Given a first principles model for the effective pairinteraction potential, U_(i,j)(r), between the two particles, thepairwise distribution function, g_(i,j), and associated function,h_(i,j), can be computed from the Boltzmann probability function in adilute suspension: $\begin{matrix}{{g_{i,j}(r)} = {{1 + {h_{i,j}(r)}} = {\exp\left\lbrack \frac{- {U_{i,j}(r)}}{k_{B}T} \right\rbrack}}} & (30)\end{matrix}$Here, the total correlation function, h_(i,j)(r) is related to the sumof all direct correlations between i and j, denoted by c_(i,j)(r) andall indirect effects on positional correlations between i and 1 mediatedthrough particle 1. One relation that relates direct and indirectcorrelations of a dense ensemble is the Ornstein-Zernike (OZ) equation(Eq 31) $\begin{matrix}{{h_{i,j}(r)} = {{c_{i,j}(r)} + {\sum\limits_{l = 1}{\rho_{l}{\int{{h_{i,l}\left( r^{\prime} \right)}{c_{l,j}\left( {{r - r^{\prime}}} \right)}{\mathbb{d}r^{\prime}}}}}}}} & (31)\end{matrix}$which requires a closure relating the total correlation, h_(i,j)(r), tothe direct correlations, c_(i,j)(r), and the pair potential forsolution. Closure relations such as the Percus-Yevick (P-Y), Roger-Young(RY), Hypernetted chain (HNC), and mean sphere approximation (MSA) mustbe employed to solve for the function h_(i,j)(r). The three dimensionalFourier transform of h_(i,j)(r) is H_(i,j)(q), where q (given by q=4πm/λ sin(θ/2)) is the magnitude of the wavevector, q. Here, m is therefractive index of the medium and θ is the scattering angle. From thesolution for h_(i,j)(r) and therefore H_(i,j)(q), the partial structurefactor function, S_(i,j)(q) can be subsequently defined as:S _(i,j)(q)=δ_(i,j) +ρ{square root}{square root over (y _(i) y _(j) )}H_(i,j)(q)   (32)where y_(i) is the number fraction of component i and δ_(i,j) is theKronecker delta function which equals unity when i=j.

In this example, we demonstrate how FDPM can be used to determineparameters which impact electrostatic interaction by considering amonodisperse suspension of known volume fraction and diameter. Theinteraction potential between charged monodisperse colloidal particlescan be expressed as a hard sphere interaction with a Yukawa tail (HSY):$\begin{matrix}{{u(r)} = \left\{ \begin{matrix}\infty & {r < \sigma} \\{{e^{2}z^{2}} < \sigma > {{\exp\left( {- {\kappa\left( {r - \sigma} \right)}} \right)}/\left( {4\quad\pi\quad ɛ\quad ɛ_{o}r} \right)}} & {r > \sigma}\end{matrix} \right.} & (33)\end{matrix}$where e is electron charge; so is the electric permittivity of vacuum;and ε is the dielectric constant of the suspending medium. The parameterσ is the particle diameter; z is the effective particle surface charge;κ is the inverse Debye screening length (which changes with the saltconcentration); and r is the center-to-center distance between twoparticles. The above HSY potential model, also termed as theone-component model (OCM), does not conserve neutrality. In it, thecharged particles are modeled as macro ions immersed in a neutralbackground but surrounded by the screened electric double layer ofthickness κ⁻¹, which decreases with increasing solvent ionic strength.The charged particles repulsively interact through volume exclusion anddouble layer overlapping. The solution of the O-Z equation using MSAapproximation with HSY potential is possible, as shown by Ginoza andYasutomi and others [Herrera, J. N.; Cummings, P. T.; Ruiz-Estrada, H.R. “Static structure factor for simple liquid metals,” Mol. Phys. 1999,5, 835-847; Ginoza, M.; Yasutomi, M. “Analytical model of the staticstructure factor of a colloidal dispersion: interaction polydispersityeffect,” Mol. P-hys. 1998, 93(3), 399-404.; Ginoza, M.; Yasutomi, M.“Analytical structure factors for colloidal fluids with size andinteraction polydispersities,” Phys. Rev. E, 1998, 58(3), 3329-3333.] Ananalytical and unique solution for the structure factor for monodispersesuspensions and the analytical partial structure factors for apolydisperse Yukawa mixture have been developed by Ginoza, M.; Yasutomi,M. “Analytical model of the static structure factor of a colloidaldispersion: interaction polydispersity effect,” Mol. Phys. 1998, 93(3),399-404.; Ginoza, M.; Yasutomi, M. “Analytical structure factors forcolloidal fluids with size and interaction polydispersities,” Phys. Rev.E, 1998, 58(3), 3329-3333.

For polystyrene latex particles, the surface charge, z, corresponds tosurface ionizable functional groups and cannot completely describe theeffective electrostatic interactions when counter ion associationoccurs. Nonetheless, the effective surface charge, z, is typically usedto compute the interaction between particles even though counter ionassociation may provide unmodeled mediation of the electrostaticinteraction. In this contribution, the effective surface charge, z, fordense colloidal suspensions is obtained by regressing experimentalvalues of isotropic scattering coefficients at varying volume fractionsto predictions made using Mie theory and Herrera's solutions for MSA-HSYinteractions.

From the FDPM measured isotropic scattering data presented in FIG. 16,the effective surface charge can be determined as illustrated in Table5. This work shows that as the salt concentration or the “screening”changes, the effective charge by which the particles experienceelectrostatic interaction changes. TABLE 5 Ionic strength, NaClequivalents 120 mM 60 mM 25 mM 5 mM 1 mM Debye screening 0.88 1.24 1.934.31 9.63 Length (nm) Effective 687 nm 1.9 1240 836 404 137 Charge 828nm 0 1005 1137 489 184 MSA-HSY

To demonstrate that the approach was sensitively evaluatingelectrostatic interactions, the salt concentration was held constant inthe dense suspension, but cationic surfactant (Rhodamine 6G), whichadsorbed onto the surface of negatively charged polystyrene, was addedin small amounts. From the scattering data, the effective surface chargedetermined from the scattering data decreased as one would expect,further demonstrating that FDPM measurements in conjunction with a firstprinciples model for particle interaction can be used to probe thenature of interactions in dense suspensions. FIG. 17 illustrates theresults of effective surface charge of the particles as a function ofcationic concentration in the dense suspension. The effective surfacecharge was determined from FDPM measurements of a 143 nm diameterpolystyrene suspension of 18.9% by volume solids with 5 mM NaClequivalents. As the concentration of cationic Rhodamine 6G surfactantincreased, the scattering increased and the effective surface charge wasobtained by regressing FDPM data for predictions using the MSA-HSY modelwith effective charge as the only variable. For systems in which theparticle size is known from another technique, FDPM measurements can beused to determine the nature of particle interactions.

It should be understood that the apparatus and methods described hereinmay also be applied to aerosols at concentrations in which multiplescattering from suspended particles may occur.

A partial list of the symbols used herein follows:

-   f(x) particle size distribution-   F scattering form factor of particles-   F_(i,j) scattering form factor for particles i and j-   k magnitude of wave vector-   K(x, λ_(i)) kernel function used in inversing non-linear integral    equation-   m refractive index of medium-   n relative refractive index of particle to medium-   S(q) static structure factor-   S_(i,j) partial static structure factor for particles i and j-   q magnitude of scattering vector-   x particle diameter-   χ² objective function used for minimization-   λ wavelength in medium-   θ scattering angle-   φ volume fraction of solids in suspension-   μ_(a) absorption coefficient-   μ_(s)′ isotropic scattering coefficient

In the foregoing written description, the invention has been describedwith reference to specific embodiments. However, after reading thisdescription one of ordinary skill in the art appreciates that variousmodifications and changes can be made without departing from the scopeof the present invention as set forth in the claims below.

1. A method for measuring an isotropic scattering coefficient forradiation having a selected wavelength and being multiply scattered in asuspension of particles, comprising: (a) placing the suspension betweena source and a detector of the radiation, the source and the detectorbeing disposed at a first selected distance apart, the radiation havingan intensity, the intensity being modulated at a first selectedfrequency of modulation; (b) measuring at least two characteristics ofthe radiation at the source and the detector, the two characteristicsbeing selected from among AC attenuation of the radiation, DCattenuation of the radiation, and phase shift of the radiation; (c)changing the selected distance apart of the source and the detector to adifferent selected distance apart and repeating step (b); and (d) fromthe selected two characteristics and a mathematical solution of thediffusion equation for propagation of light through a multiplescattering medium expressed as a linear function of the selecteddistance, calculating the scattering coefficient.
 2. The method of claim1 further comprising, after step (d), repeating step (c) so as to changethe selected distance apart to a plurality of selected distances apartand repeating step (d) for each selected distance apart.
 3. The methodof claim 1 further comprising, after step (a), changing the selectedfrequency of modulation at least once and repeating steps (b) through(d).
 4. The method of claim 1 further comprising, after step (c),calculating a criteria parameter, Q, from the results of step (b), where${\frac{{\ln^{2}\left( {{rAC}_{rel}/r_{0}} \right)} - {\ln^{2}\left( {{rPS}_{rel}/r_{0}} \right)}}{\ln^{2}\left( {{rDC}_{rel}/r_{0}} \right)} = Q},$where AC_(rel) is the ratio of AC measured at r to that measured atr_(o), PS_(rel) is the difference in measured phase at r and r_(o), andDC_(rel) is the ratio of DC measured at r to that measured at r_(o). 5.The method of claim 1 wherein in step (c) the selected distance ischanged by moving the source or the detector so as to change theselected distance apart.
 6. The method of claim 1 wherein in step (c)the selected distance is changed by substituting an alternate source ordetector.
 7. The method of claim 2 wherein the plurality of selecteddistances is at least three selected distances.
 8. The method of claim 3wherein the step of changing the selected frequency of modulation isrepeated at least twice.
 9. The method of claim 7 further comprising thestep of performing a linear regression of a logarithm of calculatedproducts of a ratio of the selected distance to an initial distance anda relative selected characteristic, at a plurality of selecteddistances.
 10. The method of claim 3 further comprising repeating steps(b), (c) and (d) for a plurality of selected distances and frequenciesand simultaneously regressing combinations of the characteristics of theradiation selected in step (b) of claim
 1. 11. The method of claim 1wherein the suspension of particles is an aerosol.
 12. The method ofclaim 1 wherein the radiation is light.
 13. A method for determining atleast one of parameters describing a non-interacting suspension, theparameters consisting of a particle size distribution and a volumefraction, comprising: measuring an experimental scattering coefficientat a selected number of wavelengths of a radiation passing through thesuspension; formulating an integral equation for a calculated scatteringcoefficient, the integral equation having an integrand including aproduct, p(x), of the particle size distribution and the volumefraction, a form factor and a structure factor; expressing p(x) as asummation of B-spline functions; calculating the form factor from Miescattering theory; calculating an inverse of the integral equation usingleast squares optimization to equate the experimental scatteringcoefficient and the calculated scattering coefficient, stabilizing theinverse using Tikhonov regularization having a regularization parameterand a non-negativity constraint on p(x) to determine a value ofseparating particle size distribution and volume fraction by anormalized distribution constraint; and using input physical data andselected upper and lower bounds of particle size, calculating volumefraction.
 14. The method of claim 13 further comprising, using thenormalized distribution constraint, calculating the particle sizedistribution.
 15. The method of claim 13 further comprising the step ofusing an L-curve to choose the regularization parameter.
 16. The methodof claim 13 wherein calculating the inverse is performed by active-setquadratic programming.
 17. The method of claim 13 wherein theregularization parameter is second order.
 18. The method of claim 13wherein the suspension is an aerosol.
 19. The method of claim 13 whereinthe radiation is light.
 20. A method for determining at least one ofparameters describing an interacting suspension, the parametersconsisting of a particle size distribution, a volume fraction and aninteraction parameter, comprising: (a) measuring an experimentalscattering coefficient at a selected number of wavelengths of radiationpassing through the suspension; (b) formulating an integral equation fora calculated scattering coefficient, the integral equation having anintegrand including a product, p(x), of the particle size distributionand the volume fraction, a form factor and a structure factor; (c)inputting an assumed or newly calculated value of particle sizedistribution and an assumed or newly calculated value of the interactionparameter; (d) expressing p(x) as a summation of B-spline functions; (e)calculating the form factor from Mie scattering theory and a particleinteraction model; (f) calculating an inverse of the integral equationusing least squares optimization to equate the experimental scatteringcoefficient and the calculated scattering coefficient, stabilizing theinverse using Tikhonov regularization having a regularization parameterand a non-negativity constraint on p(x) to determine a value of (g)separating particle size distribution and volume fraction by anormalized distribution constraint; (h) using input physical data andselected upper and lower bounds of particle size, calculating an updatedvalue of volume fraction; (i) comparing a difference between the assumedor a newly calculated value of volume fraction with the updated value ofvolume fraction to determine if the difference is less than a firstselected convergence error; (j) repeating steps (d) through (i) untilthe difference is less than the selected convergence error to determinean experimental value of volume fraction and, from p(x), an experimentalvalue of the particle size distribution; (k) calculating an updatedexperimental value of the interaction parameter using the selectedinteraction model; (l) comparing a difference between the assumed or anewly calculated value of interaction parameter with the updatedcalculated value of interaction parameter to determine if the differenceis less than a second convergence error; and (m) repeating steps (d)through (k) until the difference is less than the convergence error todetermine an experimental value of the interaction parameter for thesuspension.
 21. The method of claim 20 wherein the interaction modelincludes the LMA approximation to the HSPY model.
 22. The method ofclaim 20 wherein the interaction model includes the DA approximation tothe HSPY model.
 23. The method of claim 20 further comprising the stepof using an L-curve to choose the regularization parameter.
 24. Themethod of claim 20 wherein calculating the inverse is performed byactive-set quadratic programming.
 25. The method of claim 20 wherein theregularization parameter is second order.
 26. A system for executing acomputer program to calculate at least one of parameters describing asuspension, the parameters consisting of a particle size distributionand a volume fraction, comprising: a memory device for storing thecomputer program thereon; a processor, the processor adapted to receivea request from a device for calculation of at least one of theparameters, receive input data for physical constants, upper and lowerbounds of particle size and an experimental isotropic scatteringcoefficient, and execute the computer program to calculate an inverse ofan integral equation for a calculated scattering coefficient by using aleast squares optimization to equate the experimental scatteringcoefficient and the calculated scattering coefficient to determinevolume fraction and particle size distribution of the suspension, basedon received input data; and forward the result of the calculation to anoutput device.
 27. The system of claim 26 wherein the suspension is anaerosol.
 28. A system for executing a computer program to calculate atleast one of parameters describing a suspension, the parametersconsisting of a particle size distribution, a volume fraction and aninteraction parameter, comprising: a memory device for storing thecomputer program thereon; a processor, the processor adapted to receivea request from a device for calculation of at least one of theparameters, receive in put data for physical constants, upper and lowerbounds of particle size distribution, an experimental isotropicscattering coefficient and assumed values of volume fraction andinteraction parameter, execute the computer program to calculate aninverse of an integral equation for a calculated scattering coefficientby using a least squares optimization to equate the experimentalscattering coefficient and the calculated scattering coefficient usingthe assumed values of volume fraction and interaction parameter todetermine initial values of volume fraction and particle sizedistribution of the suspension, based on the received input data, andcompare assumed and initial values of volume fraction and interactionparameter to check for a convergence and repeat the calculation untilthe convergence is obtained; and forward the result of the calculationto an output device.
 29. Apparatus for measuring volume fraction,particle size distribution or interaction of particles in a suspension,comprising: a source and a detector of a radiation adapted for locationin a sample of the suspension at a selected distance apart; electronicapparatus for supplying a modulated radiation to the source at aselected modulated frequency; and a system for receiving signals fromthe modulated radiation, executing a computer program and displayingresults of the execution.
 30. The system of claim 29 further comprisinga plurality of sources disposed at a plurality of distances from thedetector.
 31. The system of claim 29 further comprising a plurality ofdetectors disposed at a plurality of distances from the source.
 32. Thesystem of claim 29 further comprising a mechanical device for changingthe selected distance apart of the source and the detector.
 33. Thesystem of claim 29 wherein the suspension is an aerosol.
 34. A methodfor measuring a parameter governing interparticle forces in asuspension, comprising: measuring an experimental scattering coefficientat a selected number of wavelengths of a radiation passing through thesuspension; formulating an integral equation for a calculated scatteringcoefficient, the integral equation having an integrand including aparticle size distribution, a volume fraction, a form factor and astructure factor; calculating the form factor from Mie scatteringtheory; formulating the structure factor as a function of a parametergoverning an interparticle force; calculating an inverse of the integralequation using least squares optimization to equate the experimentalscattering coefficient and the calculated scattering coefficient; andcalculating the structure factor and thereby the parameter governing theinterparticle force.
 35. The method of claim 34 wherein theinterparticle force is electrostatic.
 36. The method of claim 34 whereinthe parameter is effective surface charge.